Load Model Data

Description of two models:

Two simple models are implemented in ACT-R, performing the gambling task depending on two distinctive mechanisms.

Model1: Declarative Model (version 2.0)

The model randomly evaluates an action GUESS-MORE or GUESS-LESS, then retrieving memories about prior feedback. If “WIN” memory is retrieved, the model chooses the evaluated action. If “LOSE”/“Neutral” memory is retrieved, the model chooses the alternative action.

Model2: RL Models.

The model makes decision relying on procedural memory, pressing “MORE”/“LESS” key. GUESS-MORE and GUESS-LESS productions are competing. In the end, a feedback is given, the model receives either +1/-1/0 reward to all previous productions, and the utility of all previous fired productions will be affected.

model1.raw <- read.csv("./MODEL120210424_100307_fnca_log.csv") #%>% rename(Trial=X)
model2.raw <- read.csv("./MODEL220210424_100307_fnca_log.csv") #%>% rename(Trial=X)

#model1 <- read.csv("MODEL120210420.csv")
#model2 <- read.csv("MODEL220210420.csv")

model1.raw <- model1.raw %>% 
  mutate(CurrentResponse = case_when(Response=="j" ~ 3,
                                     Response=="f" ~ 2),
         RT = as.numeric(RT),
         Epoch = as.numeric(Epoch)) 

model2.raw <- model2.raw %>%
  mutate(CurrentResponse = case_when(Response=="j" ~ 3,
                                     Response=="f" ~ 2),
         RT = as.numeric(RT),
         Epoch = as.numeric(Epoch)) 

# only 
#model1 <- model1.raw %>% filter(ans==.1 & bll==.1 & lf==.1)
#model2 <- model2.raw %>% filter(alpha==.1 & egs==.1 & r==.1)

Compute win-stay probabilities

Win-stay probability(PSwitch) is calculated by: 1) count the number of ResponseSwitch grouped by BlockType and TrialType; 2) then calculate the proportion of ResponseSwitch (PSwitch) under groups.

future_moves <- function(responses) {
  c(responses[2:length(responses)], NA)
}
past_moves <- function(responses) {
  c(NA, responses[1:length(responses)-1])
}

count_responses <- function(model.dat) {
  model.dat <- model.dat %>%
   mutate(FutureResponse = future_moves(CurrentResponse),
         PastResponse = past_moves(CurrentResponse),
         PreviousFeedback = past_moves(TrialType),
         ResponseSwitch = if_else(CurrentResponse == FutureResponse, 0, 1))
  return(model.dat)
}

# NA for all first trial in each block
clean_previous_future <- function(model.dat) {
  model.dat <- model.dat %>%
    mutate(PastResponse=ifelse(BlockTrial==0, NA, PastResponse),
           PreviousFeedback=ifelse(BlockTrial==0, NA, PreviousFeedback),
           FutureResponse=ifelse(BlockTrial==7, NA, FutureResponse),
           ResponseSwitch=ifelse(BlockTrial==7, NA, ResponseSwitch),
           )
  return(model.dat)
}

#model1.count <- clean_previous_future(count_responses(model1.raw))

Aggregate PSwitch, RT by ParamID, TrialType, BlockType

### This aggregate looks at overall, what is the average pattern of 50 simulation (50 epoch) across different parameter set

# aggregate data by parameter set, BlockType, TrialType
m1.aggregate.ByParam.CurrentResponse.BlockType <- function(m1.dat) {
  res <- m1.dat %>%
    filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
    group_by(ans, bll, lf, BlockType, TrialType) %>%
    dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT)) %>%
    ungroup() %>%
    mutate(ParamID = group_indices(., ans, bll, lf))
  return(res)
}

m2.aggregate.ByParam.CurrentResponse.BlockType <- function(m2.dat) {
  res <- m2.dat %>%
    filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
    group_by(egs, alpha, r, BlockType, TrialType) %>%
    dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT)) %>%
    ungroup() %>%
    mutate(ParamID = group_indices(., egs, alpha, r))
  return(res)
}

# aggregate data by parameter set, BlockType, PreviousFeedback
m1.aggregate.ByParam.PreviousFeedback.BlockType <- function(m1.dat) {
  res <- m1.dat %>%
    filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback) & !is.na(CurrentResponse)) %>%
    group_by(ans, bll, lf, BlockType, PreviousFeedback) %>%
    dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))%>%
    ungroup() %>%
    mutate(ParamID = group_indices(., ans, bll, lf))
  return(res)
}

m2.aggregate.ByParam.PreviousFeedback.BlockType <- function(m2.dat) {
  res <- m2.dat %>%
    filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback) & !is.na(CurrentResponse)) %>%
    group_by(egs, alpha, r, BlockType, PreviousFeedback) %>%
    dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))%>%
    ungroup() %>%
    mutate(ParamID = group_indices(., egs, alpha, r))
  return(res)
}


#model1.aggregate <- m1.aggregate.ByParam.CurrentResponse.BlockType(model1.count)

Aggregate PSwitch, RT by Epoch, TrialType, BlockType. (This only fit for one param set)

### This aggregate looks at for single parameter set, how consistent 50 runs are

# aggregate data by Epoch, BlockType, PreviousFeedback
aggregate.ByEpoch.CurrentResponse.BlockType <- function(m.dat) {
  res <- m.dat %>%
  filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
  group_by(Epoch, BlockType, TrialType) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
  return (res)
}

aggregate.ByEpoch.PreviousFeedback.BlockType <- function(m.dat) {
  res <- m.dat %>%
  filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback) & !is.na(CurrentResponse)) %>%
  group_by(Epoch, BlockType, PreviousFeedback) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))
  return (res)
}
#model1.aggregate <- m1.aggregate.ByParam.CurrentResponse.BlockType(model1.count)

Plot functions

# plot subject aggregated data, PSwitch by TrialType x BlockType
subj.PSwitch.aggplot <- function(subj.dat, subjID) {
  res <- ggplot(subj.dat %>% filter(HCPID==subjID), 
                  aes(x = TrialType, y = PSwitch)) +
    facet_grid( ~ BlockType) +
    geom_line(aes(group = BlockType), alpha=.5) + 
    geom_point(alpha=.5, size=3, aes(col=TrialType)) +
    ylim(0,1) +
    theme_pander()
  return (res)
}

# plot model aggregated data, PSwitch by TrialType x BlockType x ParamID/Epoch
m.PSwitch.aggplot <- function(m.dat, byParam) {
  if (byParam=="ParamID") {
    res <- ggplot(m.dat, aes(x = TrialType, y = PSwitch, col = ParamID)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = ParamID), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } 
  if (byParam=="Epoch") {
    res <- ggplot(m.dat, aes(x = TrialType, y = PSwitch, col = Epoch)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = Epoch), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  }
  return (res)
}

# plot model aggregated data, T by TrialType x BlockType x ParamID/Epoch
m.RT.aggplot <- function(m.dat, byParam) {
  if (byParam=="ParamID") {
    res <- ggplot(m.dat, aes(x = TrialType, y = RT, col = ParamID)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = ParamID), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,3500) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } 
  if (byParam=="Epoch") {
    res <- ggplot(m.dat, aes(x = TrialType, y = RT, col = Epoch)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = Epoch), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,3500) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } 
  return (res)
}


#m.PSwitch.aggplot(model1.aggregate, "ParamID")
#m.RT.aggplot(model1.aggregate, "ParamID")
model1.count <- clean_previous_future(count_responses(model1.raw))
model2.count <- clean_previous_future(count_responses(model2.raw))


model1.aggregate <- m1.aggregate.ByParam.CurrentResponse.BlockType(model1.count)
model2.aggregate <- m2.aggregate.ByParam.CurrentResponse.BlockType(model2.count)

The following plots show the mean PSwitch and mean RT as a function of TrialType, BlockType. Each blue line represents a parameter set.

m.PSwitch.aggplot(model1.aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"

m.PSwitch.aggplot(model2.aggregate, "ParamID") + ggtitle("MODEL2")
## [1] "Choose one: ParamID or Epoch"

m.RT.aggplot(model1.aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"

m.RT.aggplot(model2.aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"

Summary of Pswitch vs. CurrentTrial

Model1:

Model2:

Summary of RT vs. CurrentTrial

Model1:

Model2:

model1.pf_aggregate <- m1.aggregate.ByParam.PreviousFeedback.BlockType(model1.count)
model2.pf_aggregate <- m2.aggregate.ByParam.PreviousFeedback.BlockType(model2.count)
# plot subject aggregated data, PSwitch by PreviousFeedback x BlockType
subj.PSwitch.pfaggplot <- function(subj.dat, subjID) {
  res <- ggplot(subj.dat %>% filter(HCPID==subjID), 
                  aes(x = PreviousFeedback, y = PSwitch)) +
    facet_grid( ~ BlockType) +
    geom_line(aes(group = BlockType), alpha=.5) + 
    geom_point(alpha=.5, size=3, aes(col=PreviousFeedback)) +
    ylim(0,1) +
    theme_pander()
  return (res)
}

# plot model aggregated data, PSwitch by PreviousFeedback x BlockType x ParamID/Epoch
m.PSwitch.pfaggplot <- function(m.dat, byParam) {
  if (byParam=="ParamID") {
    res <- ggplot(m.dat, aes(x = PreviousFeedback, y = PSwitch, col = ParamID)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = ParamID), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } 
  if (byParam=="Epoch") {
    res <- ggplot(model1.aggregate, aes(x = PreviousFeedback, y = PSwitch, col = Epoch)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = Epoch), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } else {
    print("Choose one: ParamID or Epoch")
  }
  return (res)
}

# plot model aggregated data, T by TrialType x BlockType x ParamID/Epoch
m.RT.pfaggplot <- function(m.dat, byParam) {
  if (byParam=="ParamID") {
    res <- ggplot(m.dat, aes(x = PreviousFeedback, y = RT, col = ParamID)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = ParamID), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,3500) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } 
  if (byParam=="Epoch") {
    res <- ggplot(model1.aggregate, aes(x = PreviousFeedback, y = RT, col = Epoch)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = Epoch), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,3500) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander()
  } else {
    print("Choose one: ParamID or Epoch")
  }
  return (res)
}

The following plots show the mean PSwitch and mean RT as a function of PreviousFeedback, BlockType. Each blue line represents a parameter set.

m.PSwitch.pfaggplot(model1.pf_aggregate, "ParamID") + ggtitle("MODEL1")
## [1] "Choose one: ParamID or Epoch"

m.PSwitch.pfaggplot(model2.pf_aggregate, "ParamID") + ggtitle("MODEL2")
## [1] "Choose one: ParamID or Epoch"

aov(PSwitch ~ BlockType * PreviousFeedback, model2.pf_aggregate  %>% filter(PreviousFeedback!="Neutral"))  %>%
  anova_summary() %>%
  xtable() %>%
  kable(digits=4) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

aov(RT ~ BlockType * PreviousFeedback, model2.pf_aggregate  %>% filter(PreviousFeedback!="Neutral"))  %>%
  anova_summary() %>%
  xtable() %>%
  kable(digits=4) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Summary of Pswicth vs. PreviousTrial

Model1:

Model2:

Summary of RT vs. PreviousTrial

Model1:

Model2:

blockId <- sort(rep(1:4, 8))

gambling <- gambling %>% 
  group_by(HCPID, RunNumber) %>%
  mutate(BlockId = blockId, 
         Response = if_else(is.na(QuestionMark.RESP), 0, QuestionMark.RESP))

bg <- gambling %>%
  group_by(BlockId, HCPID, RunNumber) %>%
  summarise(BlockType = paste("Mostly", Mode(TrialType)))

gambling <- inner_join(gambling, bg)

past_moves <- function(responses) {
  c(NA, responses[1:length(responses)-1])
}

gambling <- gambling %>%
  group_by(RunNumber, BlockId) %>%
  mutate(CurrentResponse = Response,
         FutureResponse = future_moves(Response),
         PastResponse = past_moves(Response),
         PreviousFeedback = past_moves(TrialType),
         RT = QuestionMark.RT/1200) %>%
  filter(FutureResponse != 0) %>%
  #filter(CurrentResponse != 0) %>%
  mutate(ResponseSwitch = if_else(CurrentResponse == FutureResponse, 0, 1))   %>%
  select(HCPID, RunNumber, BlockId, RT, Response, PastResponse,
         CurrentResponse, FutureResponse, ResponseSwitch, BlockType, TrialType,
         PreviousFeedback, FeedbackImage)

gambling$PreviousFeedback <- as_factor(gambling$PreviousFeedback)
gambling$PreviousFeedback <- revalue(gambling$PreviousFeedback,
                                         c("1"="Reward", 
                                           "2"="Punishment", 
                                           "3"="Neutral"))
subj1.aggregate <- aggregate.CurrentTrial(gambling)
subj1.pf_aggregate <- aggregate.PreviousTrial(gambling)

PSwitch.aggplot(subj1.aggregate) + ggtitle("SUBJ - 100307")
RT.aggplot(subj1.aggregate) + ggtitle("SUBJ - 100307")

PSwitch.pfaggplot(subj1.pf_aggregate) + ggtitle("SUBJ - 100307")
RT.pfaggplot(subj1.pf_aggregate) + ggtitle("SUBJ - 100307")
gambling <- read_tsv("../../gambling_data.txt")
gambling.trials <- gambling %>% 
  #filter(RunNumber==1) %>%
  select(HCPID, Trial, RunNumber, TrialType, RunTrialNumber, Block, 
         QuestionMark.RESP, QuestionMark.ACC, QuestionMark.RT,
         R1MostlyReward1, R1MostlyReward2, R1MostlyPunishment1, R1MostlyPunishment3, QuestionMark.OnsetDelay, 
         ConsecSameResp, ConsecLargerGuesses, ConsecSmallerGuesses, ConsecRTLess200) %>% 
  rename(CurrentResponse=QuestionMark.RESP, RT=QuestionMark.RT)  %>%
  mutate(FutureResponse = future_moves(CurrentResponse),
         PastResponse = past_moves(CurrentResponse), 
         ResponseSwitch = if_else(CurrentResponse == FutureResponse, 0, 1),
         PreviousFeedback = past_moves(as.character(TrialType)),
         BlockType = case_when(
           !is.na(R1MostlyReward1) ~ "MostlyReward", 
           !is.na(R1MostlyReward2) ~ "MostlyReward", 
           !is.na(R1MostlyPunishment1) ~ "MostlyPunishment", 
           !is.na(R1MostlyPunishment3) ~ "MostlyPunishment"
         )) %>%
  select(-R1MostlyReward1, -R1MostlyReward2, -R1MostlyPunishment1, -R1MostlyPunishment3) %>%
  filter(!is.na(TrialType))

#write.csv(gambling.trials, "../../gambling_clean_data.csv")
gambling.aggregate <- gambling.trials.session1 %>% 
  filter(!is.na(ResponseSwitch) & !is.na(RT)) %>%
  group_by(HCPID, TrialType, BlockType, PreviousFeedback) %>%
  summarise(PSwitch = mean(ResponseSwitch),RT = mean(RT))

ggplot(gambling.aggregate,
       aes(x = TrialType, y = PSwitch, col = TrialType)) +
  facet_grid( ~ BlockType, labeller = label_both) +
  geom_point(position = position_jitter(width = 0.1, height = 0.05),
             alpha = 0.1) +
  scale_color_brewer(palette = "Set2") +
  stat_summary(fun.data = "mean_cl_boot", col="black") +
  theme_pander()

gambling.trials <- gambling.trials.session1 %>% select(HCPID, Trial, RunTrialNumber, TrialType, BlockType)

gambling.trials.session1 %>%
  gghistogram("RT", add = "median")

Load Subject Data

gambling.clean <- read_csv("../../bin/gambling_clean_data.csv") %>% 
  select(-BlockType) %>% rename(BlockType=BlockTypeCoded)


gambling.cfagg <- gambling.clean %>%
  na.omit() %>%
  group_by(HCPID, BlockType, TrialType) %>%
  summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))

gambling.pfagg <- gambling.clean %>%
  na.omit() %>%
  group_by(HCPID, BlockType, PreviousFeedback) %>%
  summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))

# long format
gambling.cfagg <- left_join(gambling.cfagg %>% select(-RT) %>%
  spread(key = TrialType, value = PSwitch, sep = "."),
  gambling.cfagg %>% select(-PSwitch) %>%
  spread(key = TrialType, value = RT, sep = "."),
  by = c("HCPID", "BlockType"), suffix=c(".PSwitch",".RT"))

gambling.pfagg <- left_join(gambling.pfagg %>% select(-RT) %>%
  spread(key = PreviousFeedback, value = PSwitch, sep = "."),
  gambling.pfagg %>% select(-PSwitch) %>%
  spread(key = PreviousFeedback, value = RT, sep = "."),
  by = c("HCPID", "BlockType"), suffix=c(".PSwitch",".RT"))

#write_csv(gambling.cfagg, "../../bin/gambling_cfagg.csv")
#write_csv(gambling.pfagg, "../../bin/gambling_pfagg.csv")
gambling.clean <- read_csv("../../bin/gambling_clean_data.csv") %>% mutate(BlockType=BlockTypeCoded)
subjID = "100307_fnca"

#  Pswitch as a function of current trial type and block type
subj.cfaggregate = gambling.clean %>%
  filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>% 
  group_by(HCPID, BlockType, TrialType) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch, rm.na=T), RT = mean(RT, rm.na=T))

#  Pswitch as a function of previous trial type and block type
subj.pfaggregate = gambling.clean %>%
  filter(!is.na(ResponseSwitch) & !is.na(PreviousFeedback)) %>% 
  group_by(HCPID, BlockType, PreviousFeedback) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch, rm.na=T), RT = mean(RT, rm.na=T))

Load cleaned subject data and select one subj data = 100307_fnca Note: here, ResponseSwitch is 1 if: CurrentResponse != FutureResponse; is 0 if CurrentResponse==FutureResponse

# test whether prop is calculate correctly
subj.cfaggregate %>% filter(HCPID==subjID)
gambling.clean %>% filter(HCPID==subjID) %>%
  filter(!is.na(ResponseSwitch)) %>%
  group_by(BlockType, TrialType , ResponseSwitch) %>%
  dplyr::summarise(n=n()) %>%
  mutate(freq = n / sum(n)) %>%
  xtable() %>% kable(digits=4) %>% kable_styling(bootstrap_options = c("striped", "hover"))

#>> yes correct

The following plots show subject Pswitch pattern as a function of current TrialType and BlockType.

subj.PSwitch.aggplot(subj.cfaggregate, subjID) + ggtitle(paste("SUB", subjID, sep = "-"))

subj.PSwitch.pfaggplot(subj.pfaggregate, subjID) + ggtitle(paste("SUB", subjID, sep = "-"))

Calculate RMSE and Find Best Fit Parameter

Distribution of model1 parameters :ans, :bll, :lf

Note: This is not grid-search parameter fitting. The parameter space was chosen by optimization function by adding null_penalty (count of null responses)

model1.full.count <- clean_previous_future(count_responses(model1.raw))

model1.full.agg <- model1.full.count %>% 
  filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
  group_by(ans, bll, lf, BlockType, TrialType) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))

model1.full.agg %>% gghistogram(x="ans", binwidth = .05) + ggtitle("MODEL1: :ans distrirbutioin")

model1.full.agg %>% gghistogram(x="bll", binwidth = .05) + ggtitle("MODEL1: :bll distrirbutioin")

model1.full.agg %>% gghistogram(x="lf", binwidth = .05) + ggtitle("MODEL1: :lf distrirbutioin")

Distribution model2 parameter :egs, :alpha, and :r

dim(model2.raw)
## [1] 339200     11
model2.full.count <- clean_previous_future(count_responses(model2.raw))

model2.full.agg <- model2.full.count %>% 
  filter(!is.na(ResponseSwitch) & !is.na(CurrentResponse)) %>%
  group_by(egs, alpha, r, BlockType, TrialType) %>%
  dplyr::summarise(PSwitch = mean(ResponseSwitch), RT = mean(RT))

model2.full.agg %>% gghistogram(x="egs", binwidth = .05) + ggtitle("MODEL1: :egs distrirbutioin")

model2.full.agg %>% gghistogram(x="alpha", binwidth = .05) + ggtitle("MODEL1: :alpha distrirbutioin")

model2.full.agg %>% gghistogram(x="r", binwidth = .05) + ggtitle("MODEL1: :r distrirbutioin")

Next, we calalculate RMSE for each parameter combination, find the best parameter with min RMSE

find_best_parameter <- function(model, model.agg, sub.agg) {
  if (model=="model1") {
    m1.param <- unique(model.agg[,c('ans','bll','lf')])
    #m1.param=m1.param[1:10,]
    m1.param$RMSE = NA
    m1.param$RMSE = as.numeric(m1.param$RMSE)
    
    for (i in 1:nrow(m1.param)) {
      m1.agg <- model.agg %>% 
        filter(ans==m1.param[[i,'ans']] & bll==m1.param[[i,'bll']] & lf==m1.param[[i,'lf']])
        rmse_value <- rmse(m1.agg$PSwitch, sub.agg$PSwitch)
      m1.param[[i, 'RMSE']] <- rmse_value
    }
    return (m1.param)
  } else {
    m2.param <- unique(model.agg[,c('egs','alpha','r')])
    #m2.param=m2.param[1:10,]
    m2.param$RMSE = NA
    m2.param$RMSE = as.numeric(m2.param$RMSE)
    
    for (i in 1:nrow(m2.param)) {
      m2.agg <- model.agg %>% 
        filter(egs==m2.param[[i,'egs']] & alpha==m2.param[[i,'alpha']] & r==m2.param[[i,'r']])
        rmse_value <- rmse(m2.agg$PSwitch, sub.agg$PSwitch)
      m2.param[[i, 'RMSE']] <- rmse_value
    }
    return (m2.param)
  }
}

m1.bestp <- find_best_parameter("model1", model1.full.agg, subj.cfaggregate)
m2.bestp <- find_best_parameter("model2", model2.full.agg, subj.cfaggregate)
m1.bestp$model = "model1"
m2.bestp$model = "model2"
m1.bestp <- m1.bestp[order(m1.bestp$RMSE),]
m2.bestp <- m2.bestp[order(m2.bestp$RMSE),]

m1.bestp %>% head() %>% xtable() %>% kable(digits=4) %>% kable_styling(bootstrap_options = c("striped", "hover"))
ans bll lf RMSE model
1.2265 0.3263 0.757 0.2738 model1
1.2257 0.1000 0.100 0.2744 model1
1.2380 0.1000 0.100 0.2753 model1
1.2265 0.3263 0.739 0.2754 model1
3.0902 0.1000 0.100 0.2758 model1
1.2265 0.3263 0.729 0.2770 model1 
m2.bestp %>% head() %>% xtable() %>% kable(digits=4) %>% kable_styling(bootstrap_options = c("striped", "hover"))
egs alpha r RMSE model
0.7295 0.1038 0.6972 0.2730 model2
0.7295 0.1038 2.5030 0.2730 model2
2.1291 0.1122 5.4124 0.2732 model2
2.1287 0.1122 7.5503 0.2738 model2
0.7295 0.1038 1.4486 0.2739 model2
0.7295 0.1044 0.1000 0.2740 model2

By looking at min RMSE:

Best parameter chosen for model1 :ans=1.2264712 :bll=0.3262827, :lf=0.757

Best parameter chosen for model2 is :egs=0.7294902 ; alpha=0.1037721, :r=0.6972462

Load optimal simulation

Having the optimal parameter set, we can run simulation 800 times again by setting optimal parameter to the model. The plot shows the pattern of 800 simulation runs (Note: only single parameter set was loaded here)

subj1 = "./MODEL220210425_100307_fnca_best800.csv"
subj2 = "./seed/MODEL220210425_100307_fnca_best800.csv"
model2.best800 <- read.csv(subj1) %>%
  mutate(CurrentResponse = case_when(Response=="j" ~ 3,
                                     Response=="f" ~ 2),
         RT = as.numeric(RT),
         Epoch = as.numeric(Epoch)) 
model2.best800.count <- clean_previous_future(count_responses(model2.best800))
model2.best800.aggregate <- aggregate.ByEpoch.CurrentResponse.BlockType(model2.best800.count)

m.PSwitch.aggplot(model2.best800.aggregate, "Epoch") + ggtitle("MODEL2: optimal parameter 800 runs", 
                                                               subtitle = paste("SUB", subjID, sep="-"))

ggarrange(m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==0), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==1), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==2), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==3), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==4), "Epoch"),
m.PSwitch.aggplot(model2.best800.aggregate %>% filter(Epoch==5), "Epoch"))

Grid Search Examination

Next, the following plots show the effect of random parameter (:ANS and :EGS) on model behaviors.

m1.gsfiles = list.files(path = "./gs_check", 
                        pattern = paste(paste("^MODEL1.*", subjID, sep=""), "*_gs.csv$", sep=""), full.names = T)
m2.gsfiles = list.files(path = "./gs_check", 
                        pattern = paste(paste("^MODEL2.*", subjID, sep=""), "*_gs.csv$", sep=""), full.names = T)

m1.gs <- read.csv(m1.gsfiles[1]) %>% 
  mutate(CurrentResponse = case_when(Response=="j" ~ 3,
                                     Response=="f" ~ 2), RT = as.numeric(RT), 
                                     Epoch = as.numeric(Epoch)) 
m2.gs <- read.csv(m2.gsfiles[1]) %>% 
  mutate(CurrentResponse = case_when(Response=="j" ~ 3,
                                     Response=="f" ~ 2), RT = as.numeric(RT), 
                                     Epoch = as.numeric(Epoch)) 
m1.gs <- m1.aggregate.ByParam.CurrentResponse.BlockType(clean_previous_future(count_responses(m1.gs))) %>%
  mutate(ans = as.factor(ans), ParamID = as.factor(ParamID))
m2.gs <- m2.aggregate.ByParam.CurrentResponse.BlockType(clean_previous_future(count_responses(m2.gs))) %>%
  mutate(egs = as.factor(egs), ParamID = as.factor(ParamID))


ggplot(m1.gs, aes(x = TrialType, y = PSwitch, col = ans)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = ans), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander() +
  ggtitle("MODEL1: :ans effect", subtitle = paste("sub", subjID, sep="-"))

ggplot(m2.gs, aes(x = TrialType, y = PSwitch, col = egs)) +
      facet_grid( ~ BlockType) +
      geom_line(aes(group = egs), alpha=.5) + 
      geom_point(alpha=.5, size=3) +
      ylim(0,1) +
      stat_summary(fun.data = "mean_cl_boot", col="black") +
      theme_pander() +
  ggtitle("MODEL2: :egs effect", subtitle = paste("sub", subjID, sep="-"))

Calculate Log-Likelihood

The following analysis is to calculate the Log-Likelihood of P(m|x) given human data; the RMSE of model data and human data; The BIC scores calculated from both RSS and LL.

model2.best800.aggregate <- model2.best800.aggregate %>%
  group_by(BlockType, TrialType) %>%
  dplyr::summarise(PSwitch.m = mean(PSwitch), PSwitch.sd = sd(PSwitch))

subj1.aggregate <- subj.cfaggregate %>% filter(HCPID == subjID) %>% rename(PSwitch.subj = PSwitch)

# bic 
bic.ll <- function(k, n, logL) {
  return(k * log(n) - 2 * logL)
}
#= n log(RSS/n) + (p + 1) log(n)
bic.rmse <- function(k, n, rss) {
  return(n * log(rss/n) - (k + 1) * log(n))
}

# calculate log-liklihood
model2.best800.LL <- left_join(model2.best800.aggregate, subj1.aggregate) %>%
  mutate(PSwitch.z = (PSwitch.m-PSwitch.subj)/(PSwitch.sd),
         PSwitch.probz = dnorm(PSwitch.z), 
         PSwitch.logprobz = log(PSwitch.probz)) %>%
  group_by(HCPID) %>%
  dplyr::summarise(PSwitch.LL = sum(PSwitch.logprobz), RMSE = rmse(PSwitch.m, PSwitch.subj), RSS = RMSE^2 * 6) 


model2.best800.LL %>% mutate(BIC.LL = bic.ll(k=3, n=6, PSwitch.LL), BIC.rmse = bic.rmse(k=3, n=6, RSS)) %>%
  xtable() %>%
  kable(digits=4) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
HCPID PSwitch.LL RMSE RSS BIC.LL BIC.rmse
100307_fnca -8.1648 0.1824 0.1996 21.7049 -27.5868